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Abstract. We consider an interaction quench in the critical spin-1/2 Heisenberg 
XXZ chain. We numerically compute the time evolution of the two-point correlation 
functions of spin operators in the thermodynamic limit and compare the results to 
predictions obtained in the framework of the Luttinger liquid approximation. We find 
that the transverse correlation function (SJ'SJ' +e ) agrees with the Luttinger model 
prediction to a surprising level of accuracy. The agreement for the longitudinal two- 
point function (SjSj +e ) is found to be much poorer. We speculate that this difference 
between transverse and longitudinal correlations has its origin in the locality properties 
of the respective spin operator with respect to the underlying fermionic modes. 
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1. Introduction 

The non-equilibrium dynamics of isolated quantum systems represents a theoretical and 
experimental challenge that raises numerous fundamental questions with applications to 
a variety of fields of modern physics. In a so-called interaction quench , a system evolves 
unitarily from an initial state, which is the ground-state of a translationally invariant 
Hamiltonian differing from the one governing the time evolution by an experimentally 
tunable interaction parameter [1, 2], Over the last few years it has become clear that 
the quench dynamics in integrable systems is very unusual. While generic systems relax 
locally to standard Gibbs distributions at effective temperatures set by the initial energy 
density [3, 4, 5], integrable models retain a detailed memory of the initial state even 
at infinite times [6, 7, 8, 9] by virtue of having infinite sets of conserved charges. It 
has been conjectured that in integrable systems expectation values of local operators 
in the stationary state can be calculated using a generalized Gibbs ensemble (GGE)[7], 
a statistical ensemble formed from all local conserved charges (the role of locality of 
the integrals of motion has been highlighted in Refs [10, 11]). In order to understand 
and describe the behaviour of the quench dynamics in interacting integrable models 
new approaches have been developed [9, 12, 13, 14, 15, 16]. Investigations of quench 
dynamics in such models have focussed on the one-dimensional Bose gas (i.e. the Lieb- 
Liniger model) [17, 18, 19, 20, 21], the XXZ spin-chain [22, 23, 24, 25, 26], the Richardson 
model [27] and sine- and sinh-Gordon field theories [28, 16]. For the Bose gas, the 
knowledge of the overlaps between the BEG state and the eigenstates of the system with 
arbitrary interaction [19, 29] allowed for a characterization of the stationary state [19] 
which turned out to be compatible with a properly constructed GGE [30] , while the full 
time evolution of observables is analytically known only in the limit of strong interaction 
[31, 32] (but it is possible to exploit integrability to study the time evolution numerically 
[9, 19, 33]). For the XXZ spin chain, the overlaps with the Bethe eigenstates are known 
only for some classes of product states [34, 35, 36, 37, 38]. Again these overlaps allowed 
for a description of the stationary state in the gapped regime [24, 25] which, surprisingly, 
turned out to disagree with the predictions of a GGE formed by the infinite number of 
ultra-local conservation laws of the XXZ chain [22, 23]. This disagreement motivated 
many investigations aimed at understanding its origin, see e.g. [39, 40, 41, 42, 43, 44], 
which culminated in the recently proposed solution [45] based on taking into account 
quasi-local conserved charges. 

In spite of this intense research effort, relatively little is known about the full time 
dependence of observables in interacting integrable models in the thermodynamic limit. 
This is in stark contrast to non-interacting theories (or models that map onto such) 
[15, 31, 32, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61], In particular, 
there are no exact results for the behaviour after quantum quenches in the gapless phase 
of the Heisenberg XXZ model (with the exception of the XX chain, corresponding to 
non-interacting fermions [62, 63]). However, such quenches have been analyzed in a 
number of numerical works [62, 64, 65, 66]. The critical phase of the Heisenberg chain is 
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of particular interest, because in equilibrium many of its properties can be described by 
Luttinger liquid theory [67, 68, 69, 70]. An obvious question is then, to what extent, if 
any, Luttinger liquid theory can be used for the analysis of quench dynamics. Given that 
the quantum quench deposits an extensive amount of energy in the system, perturbations 
to the Luttinger liquid description as well as cutoff effects are expected to play a role. 
Before considering such issues one first ought to investigate whether the simple results 
obtained for quantum quench in the Luttinger model [46] can effectively describe the 
non-equilibrium dynamics of the spin chain, at least in some parameter regimes and/or 
specific time windows. Indeed this question has been addressed in Refs [71, 65] for the 
so-called Z-factor and in Ref. [66] for the stationary values of correlation functions in 
momentum space. 

Effects of perturbations away from the Luttinger model have been analyzed by 
renormalization group methods [72, 73]. This gives some insight in how perturbations 
lead to eventual thermalization on very long time scales. We mention that many other 
aspects of the quench dynamics of the Luttinger model have been considered in the 
literature [74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89]. The aim of our 
work is to quantify in some detail to what extent “simple” Luttinger liquid theory can 
be used to describe the time evolution of local observables after an interaction quench 
in the critical phase of the spin-1/2 Heisenberg XXZ chain. To that end we compare 
predictions of Luttinger liquid theory for spin-spin correlation functions to numerical 
results. 

1.1. The model and the quench protocol 

We consider is the one-dimensional spin-1/2 XXZ chain with Hamiltonian 

L 

77(A) = J X (SjSj+1 + SjSj +l + A SjSf +1 ), (1) 

3 = 1 

where Sf are spin operators at site j of a one dimensional chain and A parametrizes 
the exchange anisotropy. Throughout this paper we will set J = 1, which amounts to 
measuring all energies in units of J. The equilibrium properties of the system in the 
thermodynamic limit are well known: for |A| < 1 the system is quantum critical, while 
for | A | > 1 the ground state is antiferromagnetically ordered and excitations acquire an 
energy gap. 

In the following we want to study the out-of-equilibrium dynamics in the gapless 
phase. Our protocol is as follows: we prepare the system in the ground state |T 0 ) of 
the XX model (Ao = 0) and at time 7 = 0 suddenly quench the anisotropy parameter A 
to a finite value in the interval (—1,1] and then time evolve unitarily with Hamiltonian 
77(A). This protocol is often referred to as interaction quench. We note that the 
generalization of our work to an arbitrary quench A 0 —y A of the anisotropy parameter 
within the gapless phase is straightforward. In the following we present a systematic 
numerical study of the quench dynamics and compare the numerical results to analytic 
expressions for correlation functions obtained via Luttinger liquid theory (LL). 
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1.2. Organization of the manuscript 

In Sec. 2 we review the bosonization of the XXZ Hamiltonian and summarize how 
to use Luttinger liquid theory to obtain results for long-distance behaviour of spin- 
spin correlation functions after an interaction quench. In Sec. 3 we present the 
results of numerical computations of the same quantities by means the time-evolving 
block decimation algorithm. The numerical results are carefully compared with the 
bosonization predictions. Finally, in Sec. 6 we summarize our results. 


2. Bosonization of the XXZ spin chain 


In this section we summarize some key results on the low energy held theory description 
of the spin-1/2 XXZ chain in zero magnetic held, details can be found in e.g. [68]. At 
low energies the spin-1/2 Heisenberg chain can be described as a free compact boson 
perturbed by irrelevant operators 



K(d x 9) 2 + hd x <j>) 2 

A 


+ H» 


( 2 ) 


Here (f> = ip + (p is a canonical Bose held, 6 = ip — <p the dual held (p and <p are 
chiral components), v is the spin velocity, K is the so-called Luttinger parameter, 
and "Hin- denotes an infinite tower of perturbing operators that are irrelevant in the 
renormalization group sense. In zero held the values of K and v = voq (we recall that 
we have set J — 1) are known exactly from the Bethe ansatz solution of the XXZ chain 
[90] 


- it VI - A 2 ir 1 

v = -—, K =-—, (3) 

2 arccos A ' 2 it — arccos A 

where ao is the lattice spacing. In the XX limit A = 0 we have v — 1 and K = 1. The 
bosonized expressions for the lattice spin operators are [68, 69, 70] 


S' 2 ~ m - ^=d x (f)(x ) + (— l) j ai sin(V^7T(/>(a;)) + ..., (4) 

J Qi r 

Sj ~ b 0 (—iy cos (a/7 t6(x)) + ibi sin 9(x)) sin (\/47 t4>(x)) + ... , (5) 

where x = jao and m = {Sf). Using these expressions one arrives to the equilibrium 
two-point functions 
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In the absence of a magnetic held (m = 0) the various non-universal constants are 
known. In two-point functions the relevant amplitudes are 
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Their values are given by [69] 
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Here the parameter rj is related to the anisotropy A by 

1 / \ \ 1 

rj — 1-arccos(A) = 

7r 2A 


2.1. Interaction quench 

Our quench protocol is to prepare our spin chain in the ground state of the XX-chain, i.e. 
the Hamiltonian (1) with A = 0, and then time evolve the system with H( A) at times 
t > 0. This corresponds to a sudden quench of the interaction strength from A = 0 to 
A G (—1,1] at time t — 0. Assuming that a projection to the low-energy description in 
terms of a Luttinger liquid is possible, this corresponds to preparing our system in the 
ground state of the free boson theory A = 0), and then time evolve it with 'H(A). 
The observables of interest are the low-energy projections (4), (5) of the spin operators. 
Two point functions of these operators have been calculated for interaction quenches 
in the Luttinger model in Refs. [46, 47], and can be used in the case of interest here. 
In order to make our discussion self-contained we summarize the main points of the 
necessary calculations in Appendix A. The final results of these calculations are 

(l/A' 2 —1)/8 
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In the correlator {Sf (t)Sf +e (t)) we retain only the staggered term because the other 
contributions turn out to be small, while for (Sj(t)Sj +e (t)) the smooth and staggered 
terms turn out to be comparable in magnitude. As first pointed out in Ref. [46] for 
the Luttinger model, the power-law decays (13), (14) are very different from their 
equilibrium analogs. The amplitudes A x ! z and B z are non-universal and will be fixed 
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A 

K 

V 

a 

P 

-0.5 

3/2 

3a/3/8 

-5/36 

5/4 

-0.2 

1.14704 

0.868468 

-0.0599861 

0.315694 

0 

1 

1 

0 

0 

0.2 

0.886377 

1.12386 

0.0682023 

-0.214336 

0.5 

3/4 

3a/3/4 

7/36 

-7/16 


Table 1 . Luttinger liquid parameters for the values of A considered here. 


below by fitting (13), (14) to numerical results. They are related to the equilibrium 
amplitudes Aq and A\ via unknown relations involving the cutoffs of pre- and post- 
quench Hamiltonians. If the quench does not require a cutoff adjustment one would 
expect that 

A x ~ A x 0 , A z ~ A z , B z ~ 1. (15) 

We will see that these relations turn out to be approximately fulfilled, with deviations 
of the order of a few percent. 

As (13) and (14) are derived in the framework of a field theory approximation 
they are applicable to the description of lattice correlators only as long as £,vt, \d — 
2vt\ S> ao = 1. In Table 1 we summarize the numerical values of the parameters 
for the post-quench interaction strengths that we will consider in the following, i.e. 
A = -0.5, -0.2, 0.2, 0.5. 

3. Quantum quenches in the spin-1/2 XXZ chain: from free fermions to the 
quantum critical phase 

We now turn to the central part of this paper, a detailed numerical study of the quench 
dynamics of the XXZ spin-chain from Ao = 0 to a final — 1 < A < 1 in the gapless phase. 
To this end we employ the infinite time-evolving block decimation (iTEBD) algorithm 
[91] to study the dynamics induced by the post quench Hamiltonian (1). The algorithm 
is based on a matrix product state (MPS) description of one-dimensional lattice models 
and works directly in the thermodynamic limit. Compared to other algorithms like 
time-dependent density matrix renormalization group [92], the iTEBD has the great 
advantage of not introducing systematic errors due to finite size effects, since it works 
directly in the thermodynamic limit. This is made possible by two main features: (i) 
invariance under translation of the Hamiltonian; (ii) the possibility to parallelize the 
local updates of the time-evolving block decimation procedure. 

3.1. Method 

The iTEBD algorithm we are using is composed of two different parts: the first obtains 
an accurate MPS description of the initial state (in our case the ground state of the XX 
Hamiltonian) and the second deals with the time evolution. 
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The ground state of the XX Hamiltonian Hxx = H{ A = 0) is found using the 
iTEBD algorithm in imaginary time. The infinite chain is prepared in a given, simple 
state |T S ), which we choose as |\P S ) = 0 ? - gZ (| t)j + | \)j)/V%- This state has a trivial 
MPS representation of bond dimension xo — 1- We then evolve |T S ) in imaginary time 
with exp(— tH X x), implemented using the second order Suzuki-Trotter decomposition 
with imaginary time-step r = 0.01. We verified that further decreasing r does not 
affect the final result within our numerical accuracy. As is well known, since the 
imaginary-time evolution operator is not unitary, the MPS loses its canonical form 
(i.e. its normalisation is not constant). Hence this form must be restored while time 
elapses and in particular before taking the expectation value of any operator. We 
control the convergence of the imaginary time algorithm by keeping track of the energy 
density E 0 and waiting for it to become stationary (with an accuracy of 10~ 16 ). We 
repeat this procedure for several values of the bond dimension up to Xo — 128. A 
check of this procedure is provided by our best estimate of the ground-state energy 
density, which is Eq = —0.31830981. This differs from the exact value E\x = — 1 /7r by 
7 ■ 1CU 8 . This is quite satisfactory as the XX model is gapless and exhibits long range 
correlations. An exact MPS representation of the ground state in this case requires an 
infinite bond dimension (in other words, the singular values involved in the MPS show 
an algebraic decay). For our purposes the description of the initial state with Xo — 128 
is sufficiently accurate. Indeed, comparing the longitudinal and transverse spin-spin 
correlation functions with the exact analytical results up to distances of £ = 50, we 
find an extremely good agreement, with relative errors smaller than 3% and 1% for 
respectively longitudinal and transverse correlators. 

Using this MPS as our initial state, we can address the real time evolution. We again 
use a second order Suzuki-Trotter decomposition of the evolution operator exp (—idtH) 
with time step dt = 0.05 (we verified that the data are not affected by the time 
discretisation). We adapt the number of states used to describe the reduced Hilbert 
space by retaining, at each time step, all Schmidt vectors corresponding to singular 
values larger than X min = 1CU 12 , up to a maximum value xmax = 1024. The latter is 
actually reached fairly quickly, reflecting the fast growth of the entanglement entropy 
under time evolution. Indeed, it is well known that the computational complexity of 
the time evolution of a quantum system using algorithms based on MPS descriptions 
is essentially set by the growth of the bipartite entanglement. As the entanglement 
increases with time, we have to enlarge the dimension x °f the reduced Hilbert space in 
order to optimally control the truncation error. In spite of our refined adaptive choice 
of x, the truncation procedure is still the main source of error in our algorithm. For this 
reason, we are able to reach a maximum time T = 20 without significant truncation error 
(the Schmidt error coming from the iterative truncation of the Hilbert space remains 
always smaller than 2- 1CU 3 ). This is also reflected on the behaviour of the entanglement 
entropy of half system, which grows linearly for all explored times as it should after a 
global quantum quench [93, 54], 
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Figure 1. (Top) |(5f5J)| vs time, for different distances £ = j — i = 8, 12, 16, ... , 48 
(from top to bottom) and interaction strengths A > 0. Full lines are the numerical 
data, dashed lines represent the best fits to (13) on the interval 0 < t < t max = 20. 
(Bottom) The pre-factor A x obtained from the best fit for different values of t max . 
The dashed line is the equilibrium value of the amplitude Ag in Eq. (9). 


4. Transverse correlation functions (SfSf) 

We now turn to the two point correlation function of the x-component of spin. We find 
that (SfSf +e ) displays a strongly alternating behaviour in space, i.e. it is dominated by 
a contribution of the form (—l) e f(d), where f(£) is a smooth function. Given that this 
structure is correctly predicted by the Luttinger model, cf. (13), we will only analyze 
the absolute value |(5'f5'j : )| (as this makes plotting the correlation function easier). 

Our numerical results for | | are shown in Figs 1 and 2 as functions of time 

for different, fixed distances £, and in Fig. 3 as functions of i for the largest time t = 20. 
The main qualitative features of this correlator are as follows, (i) Apart from a brief, 
non universal transient, |(S'JS'J +£ |) is surprisingly well described by the Luttinger liquid 
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Figure 2. The same as in Figure 1 for A < 0. 


prediction (13). 

(ii) For all quenches the correlator exhibits a clear light-cone effect. At “short” 
times 1 -C 2 vt <C £, the correlations decay in a power law fashion ~ t~ a with exponent 
a = (1/A' 2 — l)/4, in agreement with (13). We note that the exponent a is related to 
properties of the Z-factor studied in Ref. [65]. The exponent a is positive (negative) 
for positive (negative) interactions A and indeed a qualitatively different behaviour is 
evident just by looking at Figs. 1 and 2. For late times 2 vt > £ the correlation function 
exhibits a power-law decay in time towards its asymptotic stationary value, which scales 
with the distance as £~ a ^ 1 ! 2 (cf. Eq. (13)). The algebraic decay as a function of distance 
at late times is shown in Fig. 3 for several quenches. Another feature visible in Fig. 3 is 
the interaction strength dependence of the velocity, at which correlations spread. This 
is reflected in the presence of “bumps” in the spatial dependence of correlation functions 
in Fig. 3 at the positions of the light cones £* ~ 2 vt. Inspection of the numerical results 
for different values of A shows that £*, and hence v, depends on A. At distances £ > £*, 
the correlators still display the power-law decay £~ 1 ^ 2 of the initial state. In particular 
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Figure 3. Correlation function {Sf Sf + f) at time t = 20, plotted as a function of 
the lattice distance £, for different interaction strengths A (different symbols). For 
comparison we also report the initial correlators (black squares). The dashed lines 
represent the asymptotic LL predictions ~ £ - “ -1 / 2 . 

for the quenches to A = —0.5 and A = —0.2 (where the propagation velocity is lower) 
these regions are visible in Fig. 3. 

Having established the main qualitative features of the transverse correlator, we now 
turn to a quantitative analysis of the crossover between the two regimes described by Eq. 
(13). To that end we perform linear fits of our numerical results to (13) as a function of 
time for fixed distances t. The parameters K and v are fixed by the quench parameter A 


A 

5A X 

SA Z 

8B Z 

-0.5 

6% 

60% 

5% 

-0.2 

2% 

10% 

5% 

0.2 

2% 

10% 

20% 

0.5 

6% 

6% 

15% 


Table 2. Approximate relative deviations of the fitting parameters from the 
equilibrium Luttinger values in Eqs. (9) and (11). The large relative error in A z (£) for 
A = —0.5 is due to the fact that equilibrium value is very close to zero (see Fig. 10). 
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Figure 4. Space-time scaling of | (Sf Sf) |£ Q+1 / 2 with t = j — i for different values of 
the post-quench interaction strengths A. After rescaling the correlators, the numerical 
results (grey lines) for different distances l collapse almost one on top of one another. 
The dashed black lines is the universal Luttinger liquid prediction. 


and reported in Table 1. We repeat the fit procedure for several time-windows [1 ,t max ] 
with t max = 12,16, 20. The fit parameter A x is found to depend only weakly on the 
choice of t max , see Figs 1 and 2). We find that the values of A x determined in this way 
are rather close to their “equilibrium” values (9), with deviations of the order of a few 
percent. In Table 2 we report estimates of the relative difference < 5A X = ( A x (£) — Aq)/Aq 
for all considered quenches, where we perform an average over the values of l which we 
consider large enough. On the other hand, fixing A x to its equilibrium value results in 
significantly poorer agreement between our numerical results and the Luttinger liquid 
prediction (13). 
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Figure 5. Fit of | (SfSf) | with t = j—i to the LL prediction taking into account cutoff 
effects through Eq. (A.12). The cutoff improves the agreement with the numerical data 
for short times, but not in the vicinity of the light cone. Left panel: t = 8,12,16, 20, 24. 
Right panel: t = 8,12,16, 20, 24, 28, 32. 


4-1. Scaling behaviour 

A different way of exhibiting the level of agreement between our numerical results and 
the LL predictions is presented in Fig. 4, which displays the space-time scaling behaviour 
of the transverse correlation function in terms of the rescaled variable 2 vt/l. The 
Luttinger liquid result (13) predicts that correlator for fixed A but different distances £ 
should collapse on top of the same scaling function, once it has been rescaled by £}C+ a . 
Of course, the usual restrictions £, vt, \£ — 2vt\ a, cf. section 2.1, continue to apply. 
Fig. 4 demonstrates that our numerical results indeed show a very nice data collapse 
for all considered quenches. This is consistent with the good stability of the fitting 
parameter A x shown in Figs 1 and 2. We note that the scaling plots in Fig. 4 are in 
good agreement with the LL prediction even rather close to the light cone: the smoothed 
peaks get sharper with increasing £ and they get closer and closer to the asymptotic LL 
results. The good agreement of our numerical results with the LL prediction is rather 
surprising, given that our fits have been performed with a single free parameter, A x , 
which moreover turns out to be very close to its equilibrium value. 

4-2. Cutoff effects 

In the previous section we concluded that Luttinger liquid theory provides a very good 
description of the quench dynamics of the XXZ Hamiltonian as long as £, vt, \£ — 2vt\ ^> 
do = 1. A natural question is whether retaining an adjustable cutoff parameter in the 
LL approach could lead to an improved description of the short-time dynamics and the 
behaviour close to the light cone. The effect of the cutoff in on the transverse spin 
correlator is given by Eq. (A. 12), where a is proportional but not necessarily equal 
to Oq. Using a as a fit parameter we obtain the results shown in Fig. 5. We find 
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Figure 6. {SfSf) as function of time for several separations & = j — i. Large 
oscillations appear when the light cone (represented by a vertical line of the same 
colour as the data) is approached. 


that, reassuringly, the amplitude A x is left unchanged by the introduction of the cutoff. 
Inspection of Fig. 5 reveals that fitting the value of the cutoff significantly improves 
the agreement with the data for short times, but not close to the light cone. The main 
remaining difference between our numerical results and the LL prediction takes the form 
of a small oscillatory contribution inside the light cone. This effect goes beyond the LL 
approximation, but is small in the transverse correlation function. This will no longer 
be the case for the longitudinal correlation function, to which we turn next. 


5. Longitudinal spin correlation functions 

The analysis of the longitudinal correlation function (S z S z +e ) is significantly more 
involved than that of (SfSf +e ). This is because the LL prediction (14) now involves 
several contributions of similar size, and because there are substantial effects not 
captured by simple Luttinger liquid theory. The latter are strongest close to the light 
cone at t* = £/2v and we now turn to their description. 


5.1. Vicinity of the light cone 


In Fig. 6 we show some typical results for quenches to positive (left panel) and negative 
(right panel) values of A. 

We fix the distance £ to be an odd integer in order to account for a strong even/odd 
effect present in the initial correlations 

(-l)'-l 


(S]S* +e ) 


t =o 


2tv 2 £ 2 


(16) 


We note that this effect is captured by the LL approximation (14), since the smooth 
and the staggered terms (proportional respectively to B z and A z ) are very close in 
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Figure 7. Space-time scaling of (SjSj +£ ) for A > 0. The numerical results (black 
and red lines) for different distances i generate two regular envelopes associated with 
even and odd £ respectively. 


magnitude. Inspection of Fig. 6 shows that as the light cone (represented by dashed 
vertical lines in the figure) is approached for a given £, large oscillations in time ensue. 
This behaviour is plainly not encoded in the Luttinger liquid approximation and can 
be understood qualitatively as follows. In the Luttinger liquid approximation the 
longitudinal spin correlation function exhibits a strong (quadratic) singularity at the 
light cone. The behaviour in this regime is determined by the high-energy modes, i.e. 
the ultraviolet part of the spectrum. However, at high energies the Luttinger liquid 
approximation is no longer expected to provide a good description of the Heisenberg 
model, as lattice effects (such as saddle points of the spinon dispersion) are important. 
On the other hand, one may hope that Luttinger liquid approximation is applicable 
sufficiently far away from the light cone. In order to investigate this possibility it is 
useful to consider the short-time (2 vt < £) and the long-time (2 vt > £) behaviour 
separately. 

5.2. Space-time scaling behaviour 

Before starting the quantitative analysis, we discuss the space-time scaling of this 
correlator. In the space-time scaling regime with £, t — y oo with their ratio fixed, the 
staggered and the smooth terms scale differently with the latter going like £~ 2 while the 
former like with /3 = K 2 — 1. Given the dependence of K on A, cf. Eq. (3), 

the staggered term is leading for A > 0 while the smooth contribution dominates for 
A < 0. We first consider the repulsive regime A > 0. In Fig. 7 we show a space-time 
scaling plot of (S'J5’J +£ )K /3+2 as a function of the scaled variable 2 vtji for several values 
of £. The most striking feature of the longitudinal correlations are their pronounced 
enhancement in the vicinity of the light cone, as compared to their initial values (in 
absolute terms the correlations are still small due to the factor of £ /3+2 ). This is quite 
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Figure 8. Space-time scaling of (S z S z + f) for A < 0. We observe fair data collapse 
at short times, and irregular envelopes at late times. 


different from behaviour seen in the transverse spin-spin correlation function. 

For a given parity of £ the numerical results are seen to exhibit a fair data collapse 
outside the light cone. The absence of symmetry around zero has its origin in the 
presence of both staggered and smooth contributions to the correlation function. In the 
vicinity of the light cone as well as in its interior the oscillations described above spoil 
any data collapse, but the different curves form nice, regular envelopes. 

In the attractive regime we consider the scaling behaviour of (S z S z +i )£ 2 . As is 
shown in Fig. 8, there again is a very strong enhancement of longitudinal correlations 
in the vicinity of the light cone. There is a good data collapse only for short times and 
large oscillations start playing an important role at an early stage of the time evolution 
(for A = —0.5 the oscillations start at 2 vt/l — 0.4). Moreover, we observe that the 
oscillations are much more irregular than for A > 0. 

5.3. Behaviour outside the light cone 2 vt < £ 

We now turn to a quantitative analysis of the longitudinal correlations outside the light 
cone, i.e. at short times t < £/2v. The time dependence of (S z S z +i ) is shown in the top 
two panels of Fig. 9 for quenches to A = 0.2, A = 0.5 and in Fig. 10 for quenches to 
A = —0.2, A = —0.5. The numerical results are compared to best fits to the Luttinger 
liquid prediction (14), which are shown as dashed lines. The agreement is seen to be 
satisfactory. The fitted values for the amplitudes A z and B z are shown in the middle 
and bottoms panels of Figs 9 and 10 respectively. 

Some comments on our fit procedure are in order. To avoid the strong light-cone 
singularity in the LL approximation at £ = 2 vt, we performed fits in the time window 
[1, (£ — 8£)/2v\ for several values of 5£. These are shown in the bottom four panels 
of Figs 9 and 10). We observe a quite satisfactory stability of the coefficients A z and 
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Figure 9. (Top) Short time regime of (S?Sj), for two different interaction strengths 
in the repulsive regime A > 0 and distances £ = j — i = 17, 21, 25, ..., 49. The 
numerical data are shown as full lines, while dashed lines represent best fits with the 
LL approximation (fitting with t max = (£ — 6)/2v). (Centre/Bottom) Amplitudes A~ 
and B z obtained from best fits with different values of t max . The dashed red lines are 
the equilibrium values of the amplitude, i.e. Ag = 1, Af = 0.0689 (for A = 0.2) and 
A\ = 0.1071 (for A = 0.5). 


B z against changing the time interval over which the fits are performed. This stability 
improves with increasing £, i.e. gets better in the regime where the LL approximation 
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Figure 10. (Top) Short time regime of (S z Sj) for two interaction strengths in 
the attractive regime A < 0 and distances i = j — i = 17, 21, 25, ..., 49. The 
numerical data are shown as full lines, while dashed lines represent best fits with the 
LL approximation (fitting with t max = (£ — 6)/2v for A = —0.2 and with t max = £/4u 
for A = —0.5). (Centre/Bottom) Amplitudes A z and B z obtained from fits with 
different t max . The dashed red lines are the equilibrium values of the amplitude, i.e. 
A z = 1, Af = 0.0357 (for A = -0.2) and A\ = 0.0179 (for A = -0.5). 

is expected to work best. We observe a small even/odd effect in the fitted values of A z 

and B z , which we attribute to subleading terms within Luttinger liquid theory. The 
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Figure 11. Numerical results for {SfSf) compared (for several distances £ = j — i) 
to the LL approximation with amplitudes fixed by fitting only short times. While the 
oscillatory behaviour is not captured by the LL approximation, it appears to correctly 
account for the decay of correlations in the accessible time window. 


fitted values of both A z and B z are approximately the same as in equilibrium. In Table 
2 we report some rough estimates of the maximal relative differences. 

5-4- LL approximation close to the light cone 

A natural question is to what extent the LL approximation, with amplitudes A z , B z 
determined above, continues to describe the correlation function close to, and inside, the 
light cone. In Fig. 11 we present results of such comparisons for two different quenches. 
The correlation function is seen to exhibit decaying oscillatory behaviour inside the light 
cone. It is evident that the LL approximation fails to capture the strong oscillations. 
However, in all cases the oscillations appear to be centered around the LL approximation 
and to decay towards the latter. 

5.5. Physics beyond the LL approximation inside the light cone 

We have seen that inside the light cone the longitudinal correlation function displays 
strong oscillatory behaviour, that is not captured by the LL approximation. We now 
turn to a parameterization of this effect. We have found that the numerical data is well 
described by the functional form 

{S]S z +e ) « B + C exp(-Dt) sin(2 vt + </>). (17) 

Here v is equal to the ground state velocity of sound and {B, C , D, <f>} are constants that 
are used to obtains best fits of the numerical results to (17). We use the time window 
\{l + 8t)/{ 2u),20] and then vary 8t in order to optimize the fit. As is shown in Fig. 12, 
the quality of these fits is quite good. 
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Figure 12. Large-time behaviour of (S'? SJ) with £ = j — i. The data (symbols) are 
compared with the best fit of Eq. (17). In the inset we report the DFT of (S?S?). 
which displays a peak close to vq = 2v. 

In order to quantify how well the oscillatory behaviour can be described in terms of 
a single frequency uq = 2v we have carried out a Fourier analysis of our numerical data. 
Discrete Fourier transforms (DFT) of (SfSj) are shown in the insets of Fig. 12. We see 
that the oscillation frequencies form a narrow band peaked around vq = 2v. Moreover, 
as the distance i — j — i increases, the the frequency distribution narrows and its centre 
approaches i/ 0 . We also performed fits of our numerical results to (17) treating v as a 
fitting parameter. As expected we found that the best value of v is then extremely close 
to the sound velocity. 

For quenches to A < 0, the accessible time window t G [0, 20] is too small to assess 
the applicability of (17). The reason is that the sound velocity v decreases rapidly 
with decreasing A, cf. Table 1. This in turn results in a small oscillation frequency in 
(17), and concomitantly few oscillations inside the light cone at times t < 20. At the 
same time our numerical computations are still limited by the growth of entanglement 
entropy, which exhibits only a weak dependence on the sign of A. This precludes us 
from exploring larger time windows. 

The fact that the oscillation frequency equals v = 2v gives us an indication of which 
processes may underlie this behaviour. According to the approach proposed in Ref. [9], 
the time evolution of the spin-spin correlation function for a quantum quench from the 
initial state |To) is given by 

(S?S? +e ) = jun Re £ e £ *s | SfSf +e |«F S ) . (18) 

\i>n) 

Here \fj n ) are eigenstates of the Hamiltonian and the corresponding energies, |$ s ) 
is a particular simultaneous eigenstates of local conservation laws of the spin-1/2 XXZ 
chain that describes the steady state, its energy, and = — ln('0|T o ). It was argued 
in Ref. [9] that only states with finite energy differences relative to contribute to 
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(18). Our numerical analysis suggests that the oscillatory behaviour inside the light 
cone is well characterized by a single frequency u 0 = 2v. In (18) this corresponds to 
energy eigenstates | if n ) with 

= ±2v. (19) 

The spectrum of the XXZ chain in the vicinity of the representative state |<!> s ) can be 
calculated [94] by a generalized thermodynamic Bethe Ansatz [95] . The most important 
excitations are particle-hole excitations with energy 

co a (k p ,h h ) = e a (k p )-e a (k h ) , (20) 

where the index a runs over the different types (“strings”) of allowed excitations at 
anisotropy A and e a (k) is their dressed energy. We conjecture (cf. Ref. [94]) that for the 
quenches considered here, the bandwidth of 6\(k) is approximately equal to that of the 
equilibrium spinon dispersion e s (k), and that the extrema of E\(k) occur at k — 0, ir/2, 7r 
(we are working at zero magnetization). For the zero temperature equilibrium spinon 
dispersion we have [90] 

|e s (7r/2)-e s (0)| = 2v. (21) 

These considerations suggest that the oscillatory behaviour inside the light cone 
originates in particle-hole excitations above the representative state connecting saddle 
points of ex (k) at k — 7t/ 2 and k = 0. These are high-energy degrees of freedom 
manifestly beyond the applicability of the Luttinger liquid approximation. 


5.6. Stationary behaviour 


Finally, we turn to the late time behaviour of the longitudinal correlations. The latest 
time accessible by our numerical computations is t — 20. In Fig. 13 we display the 
dependence of \{SjSj +l )\ on the separation £ at that time. As the correlator displays 
even/odd effects due to the presence of both smooth and staggered contributions we 
show the results for even and odd values of £ in separate graphs. The Luttinger liquid 
approximation (14) predicts that at late times 

£~ 2 for A < 0 , 

( 22 ) 


lim \(S*S* + f)\ 

t—>oo J J 




- 0-2 


for A > 0. 


where the exponent f3 is given in Table 1. Although our numerical results at t = 20 
exhibit small oscillations (probably non-stationary), they agree well with the asymptotic 
LL prediction for all values of A except for A = —0.5 and odd t. This disagreement is 
not surprising, because for this value of A the velocity v is small and the correlations 
at t = 20 are therefore not yet stationary for the largest shown values of l. 


6. Conclusions 

We have considered interaction quenches from the ground state of the XX model to 
the spin-1/2 Heisenberg XXZ chain (1) in the quantum critical regime. In equilibrium 
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Figure 13. Longitudinal correlation function (SfSij) at time t = 20 for several 
different interaction strengths A as a function of the distance t = j — i, which is 
taken to be even in the left panel and odd in the right one. The dashed lines are the 
asymptotic LL predictions (22). In the right panel the initial (t = 0) correlations are 
shown as back squares. The initial correlations for even £ are equal to zero. 


the low-energy physics of the post-quench Hamiltonian is described by Luttinger liquid 
theory, and we have investigated to that extent this description applies to dynamics out 
of equilibrium. We have computed two-point correlation functions of quantum spins 
in the XXZ chain and compared our results to predictions obtained by means of a 
Luttinger liquid approximation. Our work extends previous studies, which focused on 
other observables [65, 66]. Our main results are as follows. We found that on the time 
scales accessible to us the transverse correlator (SjSf +i ) is well described by LL theory. 
Small deviations are seen in the vicinity of the light cone, i.e. 2 vt ~ l. The good 
agreement between the LL approximation and the numerical results is surprising and 
surpasses our most optimistic expectations. Conversely, LL theory fails to give a good 
description of the longitudinal two-point function (SjS z +e ) except at rather short times. 
This implies that the longitudinal correlations are more strongly affected by corrections 
to LL theory. 

It is tempting to speculate that this difference in the applicability of the LL 
approximation could be related to the locality properties of the observables S :, -' z with 
respect to the elementary excitations that spread the correlations. While Sf is local in 
this sense, Sf is not. This can be seen by noting that the elementary excitations are 
dressed versions of the Jordan-Wigner fermions, and the expression of Sf in terms of the 
latter involves a string operator. It has been previously noted [56, 51] that local and non¬ 
local operators exhibit qualitatively different quench dynamics in particular examples. 
If our observations are indeed related to locality properties of observables, one would 
expect analogous behaviour in other models, whose low-energy equilibrium properties 
are described by Luttinger liquids. One such example is that of interaction quenches 
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in the Lieb-Liniger model. Here the field correlator (one particle density matrix) would 
then be well described by the LL approximation, while the density-density correlation 
would not. This is a opportune conjecture because in the Lieb-Liniger model interaction 
quenches (between two finite values of the interaction) are experimentally feasible, but 
there is still no analytic or numerical method to tackle its non-equilibrium dynamics. 

Finally, our results are relevant to the physics of prethermalization [96, 97, 98, 99, 
100, 101, 102, 103] in models with weak integrability breaking. A prethermalized regime 
has been observed in experiments [97] and described using Luttinger liquid methods. 
Our work shows that on intermediate time scales genuine high-energy effects need to 
be taken into account as well as “irrelevant” perturbations to the Luttinger liquid 
Hamiltonian [73]. 
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Appendix A. Interaction quench in a Luttinger liquid 


Our starting point is the Luttinger liquid Hamiltonian 


«(A) = 5 


dx 


K{d x 9) 2 + ~(d x <t>) 2 

A 


After carrying out a canonical transformation 


= <P + V = 


Vk 


, 9 = (p — ip = 9\/~K : 


(A.l) 


(A.2) 


the time evolution (at t > 0) of the chiral fields is described by the mode expansion 

~p ne ik n {x^vt) + h c 


p _ ^ 

<p(t, x ) = Q + wy( vt ~ x ) + ~7t= 

2L ^ v 4™ 


~ p _ ^ r 

<p(t, x) = Q + —(vt + x) + V -7= P-ne 

2A QAnn L 

n> 0 v 


—ikn(x+vt) 


+ ll.C. 


(A.3) 


Here we consider a ring of circumference L, k n = 2nn/L, and the various mode operators 
fulfil commutation relations 

{Q, p]=i = [Q, p] , \fim &] = s n , m . (A.4) 

At t = 0 our system is prepared in the ground state of (A.l) with K = 1 and 
v = Ja 0 . This state is conveniently constructed through a mode expansion of "H(A = 0), 
which reads 

r> _ i 

Ak n X 


p _ 2 

ipft = 0, x) = Q - —X + -H= 

Cp(t = 0, x) = Q + ^rX + 

2L n>0 V 4 


Trn 


[Pne iKx + h.c.] , 
p n e~ ik - x + h.c.] . 


(A.5) 
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The ground state at t — 0 is then defined through the property 


A>|o) 0 = o = P\0)„ = P|0}„. 


(A.6) 


Next, we use the fact that the mode operators in (A.5) and (A.3) are related through 
the transformation (A.2), which implies that 


/3n- 2 


_Vk 


+ Vk 


^ n+ 2 


_y/K 


- \[K 


F-r, 


(A.7) 


The zero mode operators are related by analogous expressions. We are now in a position 
to work out two point functions of local operators after our interaction quench. Let us 
consider 


Fit, x) = o<01 exp ex P ^°)) l°)° ■ ( A - 8 ) 

This can be calculated using the mode expansions (A.3), and then expressing the mode 
operators using (A.7). The resulting expression is then evaluated using 

o(0|e 7 ^" +<56 "|0)o = e 0 " 5 / 2 , (A.9) 


which is a simple consequence of (A.6). As we are ultimately interested in the 
thermodynamic limit we ignore the contributions arising from the zero modes. We 
find 


In (F{t, x)) 


-E 


n=1 


ik n x ik n x\ 


(2 - e lknX - e 


1 + K 2 
n L 8/1 2 

1 — K 2 

_j_ ( c ik n (x+2vt) e ik n (x-2vt) _ 2 e ^k n vt _j_ c ^ 


(A.10) 


In order to proceed, we regulate the remaining sum by multiplying the summand by a 
factor e~ 2rr£n . For e > 0 the sum then exists and is evaluated using 


ln(l-*) = -£;- 


n= 1 


(A.ll) 


After a short calculation we arrive at 


F(t,x) = 


r t 2 i 

1 4- 

1+K Z 

SK 2 

[(x + 2 vt) 2 + a 2 ][{x — 2 vt) 2 + a 2 } 

1 o 
a . 


[(2 vt) 2 + a 2 } 2 


16JT 2 


, (A.12) 


where a = eL is a short distance cutoff. In order to make contact with the XXZ 
spin chain, we need to set a proportional to the lattice spacing cio (because the 
cutoffs employed in the initial and final Hamiltonian are not the same), use that 
\x\, \x ± vt\, vt a 0 , and finally employ (5). This gives the result (13) quoted in 
the main text. All other correlators are evaluated in an analogous fashion. 
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